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ABSTRACT 

We present a 7 yr timing study of the 2.5 ms X-ray pulsar SAX J1808.4— 3658, an X-ray transient 
with a recurrence time of «2 yr, using data from the Rossi X-ray Timing Explorer covering 4 transient 
outbursts (1998-2005). We verify that the 401 Hz pulsation traces the spin frequency fundamental 
and not a harmonic. Substantial pulse shape variability, both stochastic and systematic, was observed 
during each outburst. Analysis of the systematic pulse shape changes suggests that, as an outburst 
dims, the X-ray "hot spot" on the pulsar surface drifts longitudinally and a second hot spot may 
appear. The overall pulse shape variability limits the ability to measure spin frequency evolution 
within a given X-ray outburst (and calls previous v measurements of this source into question), 
with typical upper limits of < 2.5 x 10"^'* Hz s~^ (2cr). However, combining data from all the 
outbursts shows with high (6 cr) significance that the pulsar is undergoing long-term spin down at a rate 
!> = (—5.6 ± 2.0) X 10~^^ Hz s^^, with most of the spin evolution occurring during X-ray quiescence. 
We discuss the possible contributions of magnetic propeller torques, magnetic dipole radiation, and 
gravitational radiation to the measured spin down, setting an upper limit oi B < 1.5 x 10® G for 
the pulsar's surface dipole magnetic field and Q/I < 5 x 10~® for the fractional mass quadrupole 
moment. We also measured an orbital period derivative of Porb = (3.5 ± 0.2) x 10~^^ s s~^. This 
surprising large Porb is reminiscent of the large and quasi-cyclic orbital period variation observed in 
the so-called "black widow" millisecond radio pulsars. This further strengthens previous speculation 
that SAX J1808.4— 3658 may turn on as a radio pulsar during quiescence. In an appendix we derive 
an improved (0.15 arcsec) source position from optical data. 

Subject headings: stars; individual (SAX J1808.4— 3658) — stars: neutron — X-rays: stars 



1. INTRODUCTION 

The growing class of accretion-powered millisecond X- 
ray pulsars discovered by the Rossi X-Ray Timing Ex- 
plorer (RXTE) has verified the hypothesis that old mil- 
lisecond pulsars obtained their rapid spins through sus- 
tained accretion in X-ray binaries. These objects pro- 
vide a versatile laboratory. The X-ray pulse shapes 
arising from the magnetically channeled accretion flow 
can constrain the compactness (and hence the equa- 
tion of state) of the neutron star. Tracking the ar- 
rival times of these X-ray pulses allows us to measure 
the pulsar spin evolution, which directly probes mag- 
netic disk accretion torque th eory in a particularly in- 
teresting regime (Psaltis & Ch akrabartvi ri999) and also 
allows exploration of torques arising from oth er mech- 
anism s such as gravitational wave emission (|Bildstenl 
|1998( ). There have been several reports of significant spin 
evolution in accreting millisecond pulsars, some with im- 
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plied torques that are dif ficult to reconcile with standard 
accre tion torque theorv (lMarkwardtll2003l : iMorgan et al.l 
[2003t iBurderi et al.l [20061 . 12007f ). However, a variety of 
effects (including limited data spans, pulse shape vari- 
ability, and non-Gaussian noise sources) can complicate 
the interpretation of these measurements. In this pa- 
per, we address these difficulties using a comprehensive 
analysis of the most extensive data set. 

Of the eight accretion-powered millisecond pulsars cur- 
rently known, the first one remains the best-studied ex- 
ample. The X-ray transient SAX J1808.4— 3658 was dis- 
covered during an out burst in 1996 by the BeppoSAX 
Wide Field Cameras (lin 't Zand et al]ll998f ). Timing 
analysis of RXTE data from a second outburst in 1998 
revealed the presence o f a 401 Hz (2.5 ms) accreting 
pulsar in a 2 hr binary (I Wiinands fc van der Kiidfigga 
IChakrabartv fc Morgai] Il998f ). The" source is a recur- 
rent X-ray transient, with subsequent ~1 month long 
X-ray outbursts detected in 2000, 2002, and 2005; it is 
the only known accreting millisecond pulsar for which 
pulsations have been detected during multiple outbursts. 
Faint quiescent X-ray emission has also been observed be- 
tween outbursts, a l though no pulsations were detected 
llStella et al.l 120001 : ICampana et~aI1 120021 : iHeinke et all 
120071) . A source distance of 3. 4-3.6 kpc is est ii nated 
from X-ray burst propert ies (|in 't Zand et al.l l2001l : 
iGall owav fc Cummingl l2006f l . The pulsar is a weakly 
mag netized neutron star ( [1-101 x 10® G at the sur- 
face; Psaltis fc Chakrabartvl[l999f ) while the mass donor 
is likely an extremely low-mass («0.05 Mq) brown dwarf 
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(|Bildsten fc Chakrabartvl IMl . SAX J1808.4-3658 
is the only source known to exhibit all three of 
the rapid X-ray variability phenomena associated with 
neutron stars in LMXBs: accretion-powered m illisec- 
ond pulsations ( Wiin ands k, van der KlisI IL998), mil- 
lisecond oscillations durin g thermonuclear X-ray bursts 
(IChakrabartv et al.ll2003f). a nd kilohertz quasi-periodic 
oscillations (jWiinands et all |2003[ ). After submitting 
this paper, we learn ed of an independent analysis by 
iDi Salvo et al.l ()2007[ ) reporting an increasing orbital pe- 
riod (see ^3J|) . 

An opti cal counterpart ha s been detected bo th during 
outburst (iRoche et a .1119981 : iGiles et al.lll999D and qui- 
escence (iHomer et al. 1200 It) . The relatively high optical 
luminosity during X-ray quiescence has led to speculation 
that the neutron st ar may be an active r adio pulsar dur- 
ing t hese intervals (jBurderi et al.l l2003t ICampana et al.l 
120041) ■ although rad io pulsations have not been detected 
( Burgav et al.ll2003l) . Transient unpulse d radio emission 
(Gaensler et al. 1999; Rupcn ct al." 2005) and an infrared 
excess (Wang et al. 2001; Grccnhill et al. 2006), both at- 
tributed to synchrotron radiation in an outflow, have 
been reported during X-ray outbursts. 

In this paper, we describe our application of 
phase-connected timing solutions for each outburst to 
study the spin history and pulse profile variability of 
SAX J1808. 4-3658, providing the first look at the evolu- 
tion of an accretion-powered X-ray pulsar. In section [21 
we outline our analysis methods, noting the difficulties 
raised by pulse profile noise and describing a new tech- 
nique to obtain a minimum- variance estimate of the spin 
phase in the presence of such noise. In section [3l we 
present the results of this analysis. In particular, we 
observe that the source is spinning down between out- 
bursts, the binary orbital period is increasing, and the 
pulse profiles change in a characteristic manner as the 
outbursts progress. Finally, in section IH we discuss the 
implications of these results to the properties of the neu- 
tron star and accretion geometry. 

2. X-RAY OBSERVATIONS AND DATA ANALYSIS 
2.1. RXTE data reduction 

The RXT E Prop ortional Counter Array (PGA; 
iJahoda et al.l 119961 ) has repeatedly observed 
SAX J1808. 4-3658, primarily during outburst. These 
observations total 307 separate pointings and an expo- 
sure time of 1,371 ks from 1998 through 2005. The PGA 
comprises five identical gas-filled proportional counter 
units (PGUs) sensitive to X-rays between 2.5 and 60 
keV. Each PGU has an effective area of 1200 cm^. It 
is uncommon for all five PGUs to be active: some are 
periodically disabled to decrease their rates of electrical 
breakdown (jJahoda et al.ll2006[) . The average number 
of active PGUs has declined as the RXTE ages, and 
most observations during the 2002 and 2005 outbursts 
of SAX J1808.4-3658 only include two or three. 

All but three of the observations of SAX J1808.4-3658 
were taken with the E_125US_64M_0_1S mode, which 
records the arrival of each photon with a time resolution 
of 122 /is and 64 energy channels covering the full range 
of the detectors. The other observations were rebinned 
to be compatible with the 122 /is resolution data; using 
higher time resolutions provides no benefit. We shifted 



the photon arrival times to the Earth's geocenter^ using 
our improved optical position of R.A. — 18'^08™27^62, 
Decl. = -36°58'43'.'3 (equinox J2000), with an uncer- 
tainty of 0'.'15. (Please refer to Appendix \K\ for details 
on this improved position.) We then applied the RXTE 
fine clock correction, which provides absolute time mea- 
surements with erro rs of less than 3.4 fxs (99% confidence; 
iJahoda et al.ll2006l ). Finally, we filtered the data to re- 
move Earth occultations, intervals of unstable pointing, 
and thermonuclear X-ray bursts. For three observations 
at the start and end of the 1998 outburst, we relaxed 
our requirement of stable pointing and included raster 
scanning data to extend our baseline for measuring the 
frequency evolution during this outburst. These obser- 
vations provided additional valid phase measurements, 
but they were not used to calculate fractional amplitudes 
since the contribution of the source and background var- 
ied as the RXTE panned across the source. Table [1] lists 
all the observations that we included in our analysis. 

We consistently used an energy cut of roughly 2-15 keV 
for our timing analysis. While the source is readily de- 
tectable in the PGA at higher energies, the background 
dominates above 15 keV, especially in the dimmer tails 
of the outbursts. Excluding these high-energy counts op- 
timized the detection of pulsations when the source was 
dim, providing a longer baseline for our timing analysis. 

2.2. Pulse timing analysis with TEMPO 

The core of this analysis resembles the work long-done 
for radio pulsars and slowly rotating X-ray pulsars. We 
first folded intervals of data according to a phase tim- 
ing model to obtain pulse profiles. We next compared 
the profile from each interval to a template profile in or- 
der to calculate the offset between the observed and the 
predicted pulse times of arrival (TO As). We then im- 
proved the initial phase model by fitting it to these TOA 
residuals. 

We used the tempo pulsar timing program®, version 
11.005, to calculate pulse arrival times from a phase 
model and to improve a phase model by fitting it to ar- 
rival time residuals. Tempo reads in a list of TO As and a 
set of parameters describing the pulsar timing model. It 
then adjusts the model to minimize the timing residuals 
between the predicted and observed arrival times. The 
output files include a revised timing model, a covariance 
matrix for the fit parameters, and a list of the timing 
residuals. Tempo also includes a predictive mode, which 
takes a timing model and generates a series of polyno- 
mial expansions that give the model's pulse arrival times 
during a specified time interval. Tempo has been a stan- 
dard tool of the radio pulsar community for decades and 
is well-tested at the microsecond-level accuracies with 
which we are measuring TOAs. 

Our timing models fit for the following parameters: the 
pulsar spin frequency and (if necessary) the first-order 
frequency derivative; the times and magnitudes of any 
instantaneous changes in the frequency; and the orbital 

^ We shifted the photons to the geocenter rather than the solar 
system barycenter since the TEMPO pulse timing program, which 
we used to fit phase models to the arrival times, was designed for 
radio timing and thus expects photon arrival times at some point 
on the Earth. Tempo itself performed the barycentric corrections 
using the quoted position. 

* |http: //www. atnf . cslro . au/research/pulsar /tempo/ 1 
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TABLE 1 

Observations analyzed for each outburst 



Data range # Time Avg. # Observation IDs 

(MJD) obs. (ks) PCUs 



1998 Apr 50914.8 - 50939.6 21 178.1 4.67 30411-01-* 

2000 Feb 51564.1 - 51601.9 38 126.8 3.74 40035-01-01-00 - 40035-01-04-01 

40035-05-02-00 - 40035-05-18-00 

2002 Oct 52562.1 - 52602.8 129 714.5 3.25 70080-01-*, 70080-02-* 

70080-03-05-00 - 70080-03-24-00 
70080-03-25-01, 70518-01-* 

2005 Jun 53523.0 - 53581.4 55 284.3 2.84 91056-01-01-01 - 91056-01-04-01 

91418-01-01-00 - 91418-01-07-00 

Note. — The ranges of observation IDs given here are for numerically sorted IDs, which 
do not always reflect temporal order. 



parameters. Our models supplied, but did not fit, the 
position of the source from Appendix |^ Because we 
fit the outbursts separately, the month of data that 
each provides was not sufficient to improve the source 
position: the position of SAX J1808.4— 3658 (in particu- 
lar, its right ascension) is degenerate with the frequency 
and frequency derivative on such timescales. 

To parametrize the orbit of SAX J1808.4-3658, we 
used tempo's ELLl binary model, which employs the 
Laplace parameters esino; and ecosw, where e is the 
eccentricity and w the longitude of periastron passage. 
This parametrization avoids the degeneracy of lo in low- 
eccentricity systems f|Deeter et al.iiil98li ) . For most of the 
fits, we held e = and solely fit the projected semima- 
jor axis OxSini, the orbital period Porb, and the time of 
ascending node^ Tg_sc- As a test of this assumption, we 
also repeated the fits allowing e to vary. It was always 
consistent with zero. 

This analysis depends critically on the accurate calcu- 
lation and processing of the TOAs. To verify our results, 
they were independently calculated using two entirely 
separate data pathways. One corrected the RXTE count 
data to the geocenter and used tempo to barycenter the 
TOAs, as described in the previous section; the other 
used the FTOOL^" faxbary to barycenter the count 
data. Independent codes were then used to divide the 
count data into 512 s intervals, fold it according to a 
phase model, and measure the pulse times of arrival. Fi- 
nally, we used both tempo and its replacement, tempo2 
(iHobbs et al.l[200l . to process the TOAs and refine the 
timing models. In all cases, the agreement between the 
final timing models was good. 

2.3. TO A calculation in the presence of profile noise 

Special care must be taken when measuring the pulse 
TOAs for rapidly rotating accretion-powered pulsars. In 
these systems, the pulse profiles exhibit variability on 
timescales of ~10 hr and longer that is well in excess 
of the Poisson noise expected from counting statistics. 
In this section, we develop a procedure to obtain a 
minimum- variance estimate of the timing residuals in the 
presence of such noise. 

We use the term "noise" with respect to spin timing 
analysis simply to mean phase variability of one or more 
harmonics that does not seem to be due to underlying 

Past pulsar timing of SAX J1808. 4-3658 uses the T90 fiducial, 
marking a time at which the mean longitude is 90° . Since its orbit 
is cir cular, Tgo = Taac + forb/4. 

|http: //heas arc. gsfc.nasa. gov/f tools/ 1 



spin frequency changes. While some of this profile vari- 
ability may in fact be quite ordered — distinctive pulse 
shape changes that occur in every outburst, for instance 
— we cannot model all of them and thus consider the 
unmodeled profile variability as "noise" from the phase- 
timing perspective. In this section, we attempt to min- 
imize the impact of such variability on the accuracy of 
our pulse arrival times by favorably weighting data from 
less-noisy harmonics; in t ^2.4l we describe a Monte Carlo 
technique to estimate its impact on the timing model 
parameters. 

To calculate the TOAs, we divided the timing data into 
512 s intervals and determined one pulse arrival time 
per interval. We chose this length because it provides 
sufficient counts to make accurate measurements in the 
dim tails of the outbursts, while it still is short enough to 
sample within the 7249 s orbital period. This is necessary 
to improve the binary model and resolve any additional 
short-timescale variability. 

For each outburst, we used tempo's predictive mode 
to generate a series of polynomial expansions predict- 
ing the times of pulse arrivals at the geocenter. These 
ephemerides are based on the revised optical posi- 
tion, our best-known orbital parameters, and a simple, 
constant-frequency spin model. Using the expansions, 
we calculated the expected phase for each photon arrival 
time. For each 512 s interval, we then divided the phases 
into n phase bins in order to create folded pulse profiles. 

We then decomposed the profiles into their Fourier 
components. For a given folded profile, let Xj desig- 
nate the number of photons in the jth phase bin, and 
A'ph — X]j=i tti^ total number of photons. The 

complex amplitude of the fcth harmonic is then 

n 

flfc = Xj exp {2TTijk/n) . (1) 

Throughout this paper, we number the harmonics such 
that the fcth harmonic is k times the frequency of the 
401 Hz fundamental. Since our analysis for the most part 
handles the phases and amplitudes of harmonics sepa- 
rately, we define these quantities explicitly as follows: 

Ak exp [2TTik {(t)k + ^A)] ^ 2ak . (2) 

Here we are interested in the amplitude^"'^ Ak and the 

We define A^. such that it is the actual amplitude, in photons, 
of the observed pulsations. We must therefore include both the 
positive and negative frequency components (which are equal for 
real signals), introducing the factor of 2 on the right-hand side of 
eqs. 1I2J and 



4 



Hartman et al. 



phase residual A(/)fc, which we measure relative to a fixed 
phase offset cfik- We include these offsets because we 
are principally interested in measuring the phase devi- 
ations from a fixed template profile, which we obtain by 
transforming the overall folded pulse profile from an out- 
burst:^^ 



A'f. exp {2nik(t)k) = 2 



2^ ^3 



exp [2'nijk/r 



(3) 



Here x'^ and N!^Yi the phase bin counts and total 
counts for the template pulse profile. 

Note that we define the phases such that shifting a 
fixed pulse profile by some phase A0 produces the same 
shift in the phase of each of its harmonic: A0fc = A(/>. 
Hence the unique phases for each A0fc range from to 
1/k. Positive phase residuals corresponding to time lags: 
A(/)fc > indicates that the fcth harmonic arrived later 
than predicted by the model. 

The uncertainty in the phase residuals A(/)fc due to 
Poisson noise (derived in Appendix [B]) are 



ph 



(4) 



For our analysis, we rejected phase measurements with 
uncertainties greater than 0.1 ms (i.e., 0.04 cycles). Gen- 
erally, this cut only removed points in the tails of the 
outbursts, where the fiux was low. 
The measured fractional rms amplitudes are 



A, 



V2{Np^-B) 



(5) 



where B is the approximate number of background events 
within our energy range and time interval, estimated us- 
ing the FTOOL pcabackest.^^ The add in quadra- 
ture: the total rms fractional amplitude for a pulse pro- 
file described with m harmonics is r = Cl2T=i''"k)^^'^ ■ 
Uncertainties on the fractional am plitudes are com- 
puted using th e meth od described bv lGrotlJ ()1975D and 
IVaughan et al.l (|1994f ). which accounts for the addition of 
noise to the complex amplitude of the signal. The prob- 
ability that the detection of a harmonic is due solely to 
Poisson noise is exp(— P/j), where Pk = ■jA'l/Nph is the 
unit-normalized power for the fcth harmonic. For a fuller 
review of Fou rier techniques in X-ray timing, we defer to 
Ivan der KlisI ([19891 . 

Since each harmonic provides an independent measure- 
ment of the phase residual, we can combine them to pro- 
vide the overall phase residual for the sample pulse. We 
obtain the optimal estimator by weighting each measure- 
ment according to its variance: 



Thus far, this analytical method closely parall els the 
work long done on spin-powered pulsars (e.g., iTavlod 
119931 Appendix A). 

However, there are some essential differences that must 
be taken into account when dealing with accretion- 
powered pulsars. Unlike spin-powered pulsars, which 
usually show one or more sharp, asymmetric pulses per 
cycle, there is little harmonic content in the pulsations 
of SAX J1808.4-3658 beyond fc = 2, so we truncate the 
series there. Additionally, while the individual pulses of 
radio pulsars show appreciable variability from one pe- 
riod to the next, their integ, rated profiles are very stable 
([Manchester fc Tavloilll977l ). In such cases, one expects 
that the pulse fractions of sample pulses are similar to 
the template (Afc/iVph « ^k/-^ph) ^^"^ that the harmon- 
ics refiect a common phase residual {A(j)k « Ac/)) that 
traces the rotational phase of the star. Indeed, the stan- 
dard template-matching analysis is predicated on these 
assumptions. Furthermore, any variability is assumed to 
be due to Poisson noise, which is of equal magnitude at 
all timescales (i.e., it is white noise). In contrast, the 
accretion-powered pulsars show substantial pulse profile 
variability. Beyond the usual Poisson noise (cfe from 
eq. ID), three additional issues complicate the usual ap- 
proach of template matching: long-timescale correlations 
(i.e., red noise) in the observed pulse fractions, with each 
harmonic's Ak varying independently; red noise in the 
phase offsets A^^; and sudden pulse profile changes, in 
which the phase offset between the two measured har- 
monics changes drastically on the timescale of the obser- 
vations. 

In their ti ming analysis of 2 83 s a ccretion-powered pul- 
sar Vela X-l. lBoynton et all ([19841 ) partially address the 
issue of intrinsic pulse profile noise in the harmonics. 
In the most general case, the amplitude of the vari- 
ability in the phase residuals Acj)^ is different for each 
harmonic k. This is the case for both Vela X-1 and 
SAX J1808.4-3658. They correct for the harmonic de- 
pendence of these fiuctuations by scaling the phase resid- 
uals of each harmonic by constants chosen such that the 
phase residuals all have t he same amplitude of variability 
([Bovnton fc Deeterlll985D . Thus the infiuence of particu- 
larly noisy harmonics was diminished, and they were able 
to measure with much greater accuracy the underlying 
spin of Vela X-1. 

Our approach was similar. For each outburst, we mea- 
sured the total rms amplitude of the phase residuals 
for each harmonic with respect to a best-fit constant- 
frequency model. These residuals will represent the com- 
bined effect of the Poisson noise and any intrinsic profile 
noise: 



(8) 



Wk ^ 



k=l 
1 



wkA4)k 



PAl 



fc=i 



Wk 



(6) 



(7) 



'k ^^Ph 

In the presence of sudden pulse profile changes, we may use 
multiple templates during a single outburst, thus using different 
values of </)j. on either side of the change. Further description is at 
the end of this section. 

|http: //heasarc.gsf c .nasa. gov/docs/xte/recipes/pcabackest .hlldiib] outbursts from skewing the results 



We calculated the Poisson contribution (a'f.) as a 
weighted mean^** of the results from equation ^ , giving 
us a value for (t| ju^. . We then incorporate this additional 

^■^ To calculate the mean of the variances cr^, we weight each 
according to equation JTJl. Labeling each uncertainty as a^i for 

intervals I = 1, . . . , A^int, we have {o'^) = (^i^i* ^^l^ / ^intj 
This weighting scheme prevents large variances during the tails of 
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uncertainty into our weighting to determine A</): 



cr^ changes from one TOA measurement to the next due 
to the variabihty of the pulse fraction and count rate; 
there is no assumption that these are constant, as there 
is in the case of standard template fitting, tr^ is a con- 
stant measured independently for each harmonic of each 
outburst. The result is a minimum- variance estimator 
for A(f>k. For instance, if the 802 Hz second harmonic 
has smaller intrinsic fluctuations than the fundamental, 
then our method presumes that it better reflects the spin 
of the NS and will weight it more strongly. 

Sudden pulse profile changes are somewhat simpler to 
deal with. We use a different template pulse profile (and 
hence different measurements on either side of the change 
of (j)k, defined in eq. [3]). We only modeled one such 
sudden profile change in this way: at the end of the 
main body of the 2002 outburst (around MJD 52576), 
the fundamental phase 0i experienced a shift while the 
second harmonic, (j)2, remained constant. The stability 
of (t>2 allowed us to phase connect across the feature, as 
iBurderi etHII (|2006f) also noted. 

Our distinction between sudden profile changes and 
pulse profile noise is admittedly somewhat arbitrary. We 
make it solely in the interest of best estimating the rota- 
tional phase of the star — we are not claiming to model 
some underlying difference in physical processes. In 2002, 
the phase residuals on either side of the modeled pulse 
profile change were quite stable, albeit with different val- 
ues of 01. This stability in both harmonics makes it 
a good candidate for such treatment. In contrast, the 
phase residuals of both harmonics during the 2005 out- 
burst show greater amplitude fluctuations at nearly all 
timescales. While its residuals and the 2002 residuals 
follow a similar pattern at the end of the main body of 
the outburst (the phases of the second harmonic remain 
roughly constant, while the phases of the fundamental 
drop appreciably) , the phase of the fundamental contin- 
ues to fluctuate wildly after this event rather than set- 
tling down on a "new" template profile. Therefore we 
elect to attribute these profile changes to intrinsic noise 
and weight the relatively stable second harmonic more 
strongly. 

However, when the phases of both harmonics are con- 
tinuously changing, the ability to define pulse arrival 
times breaks down, and the data are of little use for de- 
termining the spin of the star. For SAX J1808.4-3658, 
the pulse profile is changing throughout the rises and 
peaks of the outbursts, so we excluded these data from 
our measurements of the spin frequency. We did in- 
clude these data when calculating the orbital parame- 
ters. Since the timescale for these pulse profile changes 
(> 10 hr) is many times the orbital period, they tend 
to average out and have little impact on these measure- 
ments. 

The resulting phase residuals give the best estimator 
for the offset between the measured and predicted pulse 
arrival times. By adding these offsets to the phases pre- 
dicted by the tempo ephemerides, we arrived at more 
accurate pulse arrival times for each interval. 

2.4. Parameter fitting and uncertainty estimation 



After measuring the pulse times of arrival, we input 
them into TEMPO to refit the timing solution. In order to 
interpret the resulting models, we must understand the 
nature of the noise in the TOAs and how it affects the 
model parameters. The harmonic weighting system de- 
scribed above makes the optimal choice to mitigate the 
phase variability due to a particularly noisy harmonic, 
but the TOA residuals are still of substantially greater 
magnitude than would be expected from Poisson noise 
alone. This leaves us with the task of estimating the fit 
uncertainties in the presence of such noise. These uncer- 
tainties are crucial to our construction of timing models, 
as they are needed to estimate the significance of fit com- 
ponents, such as frequency derivatives and instantaneous 
frequency changes. 

We noted in the previous section that we treat the 
TOA residuals as noise, despite some of their variabilty 
arising from pulse shape changes that recur in every out- 
burst. This is not a bad approximation: because we fit 
our models separately for each outburst, correlations in 
the pulse profile variability between outbursts are not 
relevant. Furthermore, the power spectra of the TOA 
residuals (see Fig.|4]in §3.3p resemble the power- law noise 
spectra typically observed in actual red noise processes, 
so treating it as such is reasonable. 

We initially used the simplest possible timing model 
when fitting the TOAs of each outburst in TEMPO: a cir- 
cular orbit and a constant frequency. (Note that we fit 
independent models for each outburst. The uncertain- 
ties were too large to phase connect between outbursts.) 
When this simple model proved insufficient to account 
for the phase residuals, we introduced a nonzero and 
instantaneous frequency changes, as needed. However, 
there is a danger of overfitting the data. It is important 
to recognize that some of the features in the residuals are 
probably pulse profile variability rather than spin evolu- 
tion. We took great care in our attempts to distinguish 
between the ;> measurements and the artifacts of intrinsic 
timing noise. 

The colored nature of the timing noise in both harmon- 
ics is the primary difficulty in the interpretation of the 
parameter fits. Tempo assumes that the TOA uncertain- 
ties one gives it are white and approximately Gaussian, 
as is the case of pure Poisson noise. As a result, it sys- 
tematically underestimates the uncertainties in the fitted 
parameters in the presence of timing noise. Red timing 
noise is particularly problematic, because it dominates 
on the long timescales on which v and measurements 
depend. 

Instead of adopting this white noise assumption of 
TEMPO, we estimated confidence intervals for i/ and !> 
using Monte Carlo simulations of the timing residuals of 
each outburst. After using tempo to obtain the best 
fit for a timing model, we calculated the power spec- 
trum P(/) of the timing residuals that tempo output. 
This spectrum is a convolution of the true noise spec- 
trum and the sampling function; most notably, there is 
excess power around 1 d, an artifact of RXTE observa- 
tions often being scheduled approximately a day apart, 
and at the RXTE orbital period of 96 min due to Earth 
occultations. We applied a low-pass filter to remove these 
peaks in an attempt to approximate the underlying noise 
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spectrum, P'{f). 

P'U) = P{f) X [(1 - A) eM-fr!) + A] . (10) 

Tc gives the time scale for the fow-pass cutoff. A gives 
the fraction of high-frequency noise to let through, repro- 
ducing the short-timescale scatter (principally but not 
entirely Poisson) that we observed within each observa- 
tion. Typical values were 3 d and 10-20%. 

We then created thousands of sets of artificial phase 
residuals with the noise properties of the filtered spec- 
trum, /"(/). To reproduce the sampling irregularities, 
we removed all points at times absent in the original data. 
The parameters of the low-pass filter were tuned such 
that the mean power spectrum of the resulting Monte 
Carlo residuals was as close as possible to the original 
power spectrum, P{f). For each set of residuals, we 
measured the frequency of the best linear fit (or, if our 
TEMPO model fit for ly, the frequency derivative of the 
best quadratic fit). The standard deviations of these 
measurements provided uncertainty estimates for the re- 
spective parameter, this time more accurately accounting 
for the noise spectrum. 

We relied solely on tempo to calculate the uncertain- 
ties in the binary orbit parameters. While the intrinsic 
pulse profile noise spectrum is colored on the timescale 
of days to weeks, the phase residuals are approximately 
white on timescales equal to and shorter than the 2 hr or- 
bital period. The amplitude of the short-timescale vari- 
ability is roughly 1.5 times what one would expect from 
counting noise alone, so we scaled our Poisson-derived 
phase uncertainties accordingly when estimating the un- 
certainties of the orbital parameter fits. This rescaling 
makes tempo's uncertainty estimates for the orbital pa- 
rameters reasonably accurate. One important consis- 
tency check of this simple approach worked nicely: the 
1998, 2002, and 2005 measurements of Porb and CxSini, 
two parameters that should be the same for each out- 
burst at our level of accuracy, were indeed found to be 
constant, with reduced statistics close to unity. 

In deriving new binary and spin parameters, the new 
values sometimes differed considerably from the param- 
eters with which we initially folded the data for TOA 
calculations. If our orbital model improved substantially, 
we iterated the above procedure, calculating new times of 
arrival for each 512 s interval and refitting. Because the 
orbit introduces a periodic frequency modulation with 
amplitude Av = vq ■ 27rax sin«/cPorb > 1/512 s, an inac- 
curate orbital ephemeris can significantly reduce detec- 
tion strength. In contrast, the spin frequency is remark- 
ably stable through all the observations, so there was 
no need to recalculate TOAs upon the relatively minor 
revisions to the spin model. 

3. RESULTS 

The results of our pulse timing solutions are shown in 
Figure [1] which compares the light curves, phase resid- 
uals, and fractional amplitudes for each outburst. In- 
specting the best-fit frequency lines in the phase residual 
plots, it is clear that a constant pulse profile attached 
to a constant-frequency rotator does not adequately de- 
scribe the observed residuals. We consider five sources of 
phase residuals relative to a best-fit constant-frequency 
model: Poisson timing noise, intrinsic pulse profile noise, 
sudden and well-defined pulse profile changes, additional 



spin frequency derivatives, and instantaneous frequency 
changes in the underlying rotation of the star. In this 
section, we will consider all these possible contributions 
to the residuals and their relationships with each other 
and the other properties of each outburst. 



3.1. Light curves of the outbursts 

The light curves of each outburst are quite similar 
in shape. We divide them into four stages: the rise, 
which was only definitively captured in 2005 and took 
«5 d; the short-lived peak at a 2-25 keV flux of (1.9- 
2.6) X 10~^ erg cm^^ s~^, equal to a luminosity of (4.7- 
6.4) X 10^^ erg s~^ using the distance of 3.5 kpc and 
bo lometric correction of Lboi/^ 2-25 kcV = 2.12 derived 
bv lGallowav fc Gumming! (|2006( l: a slow decay^^ in lumi- 
nosity, lasting 10-15 d, until the source reaches approxi- 
mately 8x 10~^° erg cm~^ s~^(= 2.0 x 10^^ erg s^^); and 
a sudden drop followed by low-luminosity flaring as the 
outburst flickers out, with the timescale between flares 
on the order of 5 d. Figure [2] shows a cartoon of a typi- 
cal outburst from SAX J1808.4— 3658, with each of these 
stages labeled. 

The RXTE first collected high-resolution timing data 
from SAX J1808.4-3658 during the 1998 Aprfl outburst. 
Data from the RXTE AU Sky Monitor (ASM) show 
that the peak luminosity occurred approximately three 
days before the first PCA observation. (See Fig. 2 of 
iGaflowav fc Ciimmind[2006l for a comparison of the ASM 
and PCA light curves; note that its 1998 plot does not 
include two raster scans analyzed here.) Unfortunately, 
PCA observations stopped shortly after the body of the 
outburst and do not sample the tail. 

SAX J1808.4— 3658 was discovered to be again in out- 
burst when it emerged from behind the Sun in 2000 Jan- 
uary. Coverage of the o utburst was limi t ed an d included 
only the outburst tail. IWiinands et al.l (j2001l ) comment 
on the erratic nature of the flaring during the tail. ASM 
data indicate that the peak occurred 2 weeks prior to 
the first PCA observation and that we are observing the 
later, dimmer stage of the flaring tail. Comparison of 
the PCA data with the 2002 and 2005 outbursts suggests 
likewise. 

The 2002 outburst, detected in mid-October and ob- 
served for the next two months, was the brightest, had 
the best PCA coverage, and included the detection of 
four extremely bright thermonuclear X-ray bursts dur- 
ing its peak. Its light curve was very similar in shape to 
the 1998 outburst. 

In 2005 June, SAX J1808.4-3658 was again in out- 
burst. This time, the detection preceded the peak by a 
few days, providing a full sampling of the light curve. 
This outburst was somewhat dimmer, with a peak lumi- 
nosity of only 70% of the 2002 peak and a correspond- 
ingly shorter slow-decay stage. The subsequent rapid 
decay and flaring tail look quite similar to the other out- 
bursts. 

15 Some authors ('e.g.. lCui eraI|[T99g: FBurderi et aI.II20061) refer 
to this part of the outburst as the "exponential decay" stage, based 
on the approximately exponential dimming of the 1998 and 2002 
outbursts. However, the fall off of the 2005 outburst during this 
stage is closer to linear, so we simply refer to it as the "slow decay" 
stage to contrast it with the more rapid luminosity drop at its 
conclusion. 
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Fig. 1. — The light curves, phase residuals, and fractional amplitudes for all four outbursts. The top panels show the background-subtracted light curves for each outburst. The strips 
along the top of the graphs indicate the times of observations; stars indicate the times of thermonuclear X-ray bursts. The second panels show show the phase residuals relative to a 
constant frequency of 400.97521025 Hz, with black points giving the phases of the fundamental and grey points indicating the second harmonic. Positive phases indicate pulse arrivals 
later than predicted by the phase model. The error bars reflect the statistical errors only, as calculated in equation (|4j. These data have been binned so that there is only one point per 
observation. The black lines indicate the best-fit constant-frequency model for each outburst. The bottom plot shows the fractional amplitudes of the fundamental and harmonic. 
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Fig. 2. — The anatomy of a typical outburst from 
SAX J1808.4-3658. The features of the hght curve and their 
fiuxes and timescales are similar to those observed during the 
1998, 2002, and 2005 outbursts. Bolometric luminosities assume 
a distance of 3.5 kpc and a correction of Lbol/i2— 25 keV = 2.12 
UGallowav fc Cummin3l200^') . 

3.2. Characteristic pulse profile changes 

Just as the light curves of each outburst were quite 
similar, the evolution of the pulse profile during each 
outburst was remarkably consistent. Figure [3] illustrates 
the full range of pulse profiles that we observed from 
SAX J1808.4— 3658. In many instances, the similarity 
of the pulse profiles between outbursts is quite striking. 
In this section, we describe how these profiles change 
throughout the outbursts. 

We observed the outburst rise, labeled as profile 1 in 
Figure [31 exclusively during the 2005 outburst. The pro- 
files are smooth and asymmetric, with a slow rise followed 
by a more rapid drop-off after the peak. There is no sign 
of a second peak. 

We observed the outburst maxima during 2002 and 
2005. The similarity of the pulse profile evolution be- 
tween the two outbursts is remarkable. During the first 
half of the maxima (labeled as profile set 2), the profiles 
show a secondary bump lagging the main pulse. Com- 
pared to the burst rise, the fractional amplitude has de- 
creased somewhat. During the second half of the maxi- 
mum, the pulse becomes broader, subsuming the lagging 
secondary bump. (See profile set 3.) This change appears 
to be gradual: in both outbursts, a mid-peak observation 
exhibited an intermediate pulse profile. 

Profile set 4 shows the pulse profiles during the slow 
decay stages of the outbursts. The 1998 and 2002 profiles 
are quite similar: the pulses are somewhat asymmetric, 
rising more steeply than they fall. In both outbursts, 
this pulse profile is very stable during the approximately 
10 d of the decay in luminosity. During the 2005 out- 
burst, this asymmetry is more pronounced, and the pro- 
file varies between observations. Initially, the pulse ex- 
hibited a small lagging bump (profile 4A), quite similar 
to the pulse profile during the first half of the outburst 
maximum. The relative size of that bump varied sub- 



stantially, in some observations appearing as a small sec- 
ondary peak (profile 4B). Over the course of the decline, 
the source switched back and forth between a double- 
peaked and single-peaked profile as indicated in the fig- 
ure. A given state would typically be seen for two or 
three observations (1-2 d) before switching to the other. 

Profile set 5 covers the rapid drop in flux at the end of 
the outbursts. During 1998, the pulse profile was quite 
stable and did not appreciably change during this drop, 
although its fractional amplitude increased somewhat. 
In contrast, the 2002 and 2005 outbursts show a major 
pulse profile shift concurrent with the drop in luminosity. 
Prior to the drop, the pulses in set 4 show a quick rise 
and a slower fall. After the drop, the asymmetry of the 
2002 and 2005 profiles reverses: profiles 5B show a slow 
rise and a quick drop. In terms of harmonic components, 
these changes represent a shift in the phase of the fun- 
damental by approximately 0.15 cycles as it went from 
leading the second harmonic to lagging behind it. The 
phase of the harmonic did not change. In both outbursts, 
observations during the ~2 d of rapid luminosity decline 
reveal an intermediate stage in which the main pulse is 
momentarily symmetric (profiles 5A). During this tran- 
sition, small but significant secondary pulses are present. 

During the flaring tail of the outburst (profile set 6), 
the pulse proflle again showed substantial variability. In 
2002, the proflle repeatedly switched between an asym- 
metric pulse (profile 6A, identical to the pulse profile at 
the end of the rapid dimming stage) and a double-peaked 
profile (profile 6B). The double peaked pulse profile oc- 
curs principally (but not exclusively) at the end of the 
fiares, as their luminosity declines. These pulse profile 
changes are almost entirely the result of changing frac- 
tional amplitudes of the harmonic components; the phase 
offset between the fundamental and second harmonic re- 
mains for the most part constant. A notable exception 
occurs during the decay of the first flare at around MJD 
52582. At this time the phase of the fundamental jumped 
by «0.2 cycles, indicating a sudden lag of this amount 
behind its previous arrival time. By the next observa- 
tion, less than two hours later, the phase residual of the 
fundamental returned to its previous value. 

The tail of the 2005 outburst is more chaotic. The 
fractional amplitudes and phases both exhibit strong red 
noise, producing a pulse proflle that is sometimes asym- 
metric with a slow rise and quick fall (6A); at other 
times asymmetric with a quick rise and slow fall (6C; 
not shown, but basically just the reverse of proflle 6A); 
and in one instance clearly double-peaked (6B). The 
observations were sparse and generally short, so it was 
impossible to better characterize the evolution of these 
pulse proflle fluctuations. The flaring tail of the 2000 
outburst was quite similar, with a highly variable pulse 
proflle that included double-peaked proflles and asym- 
metric single pulses of both orientations. We did not 
include it because the observations were few and sparse. 

3.3. Noise properties of the timing residuals 

To measure the spin phase of SAX J1808.4— 3658 using 
the formalism developed in ^2.'i[ we must characterize the 
variability of the harmonic components. This variability 
encompasses both the pulse proflle changes discussed in 
the previous section as well as any noise in the spin phase 
of the star. 
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Fig. 3. — A comprehensive view of the 2-15 keV pulse profiles observed from SAX J1808.4— 3658. Each pulse profile was calculated by 
folding the observations within the indicated time intervals using the best-fit constant-frequency model of each outburst, so any movement 
of the peaks reflects the phase offsets from the constant frequency. The profiles are background-subtracted, normalized such that the phase 
bins have a mean value of unity, and plotted on 0.80-1.20. Thus the plotted profiles accurately show the change in fractional amplitude 
during the outburst. The profiles are numbered according their position within the outburst: 1 indicates the burst rise; 2, the beginning 
of the outburst maximum; 3, the end of the maximum; 4, the slow decay stage; 5, the steep luminosity drop marking the end of the main 
outburst; and 6, the flaring tail. During some parts of the burst, two pulse profiles are present, with the source switching between them. 
In these cases, we show both profiles and label the regions of the light curve in which they occurred accordingly. The solid black line shows 
the fluxes from the PCA observations; the grey boxes show the fluxes from the ASM daily averages. 



In our analysis of the phase residuals, we took into 
account the rms amplitude of the intrinsic pulse profile 



noise in each harmonic, a 



fe.iiit ' 



defined in equation 



A casual glance at the phase residuals of Figure [T] re- 
veals that the magnitudes of these fluctuations vary sub- 
stantially between outbursts. Table [2] summarizes these 
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Fig. 4. — Power spectra of the phase residuals for the fundamen- 
tal (black squares) and the second harmonic (grey circles) relative 
to the best-fit constant-frequency models. (The 2002 model also 
includes a phase shift in the fundamental to account for the profile 
change at MJD 52577.) The dashed lines show the power level 
due to counting statistics, a white-noise contribution proportional 
to (o"^). The data points show the powers Pk{f), from which 
we nave subtracted the contribution of counting statistics. These 

powers are normalized such that J^q-y Pk{f) df = (""fc int)' 
defined in equation Jsjl . The vertical dotted lines show the relevant 
time scales for the spectra: the 96 min and ftil d periodicities of the 
RXTE observations, and the 121 min SAX J1808.4-3658 orbital 
period. 



amplitudes for each outburst and compares them to the 
mean amphtudes of their Poisson noise, (cr^)^^^- These 
values are then used in equation ([9]). For instance, in the 
2002 outburst the fundamental is more heavily weighted 
in measuring the spin phase than the second harmonic, 
while in 2005 the opposite is true. 

The scatter of the phase residuals between observa- 
tions is generally greater than the scatter within an ob- 
servation, suggesting that the pulse profile noise is red. 
Power spectra of the phase residuals, shown in Figure [H 
confirm this. We estimated these power spectra using 
Fourier transforms of the residuals from equally spaced 
512 s bins. Here we have not attempted to deconvolve 
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Fig. 5. — The fractional amplitude of the second harmonic scales 
with flux according to a power law of slope —0.50 ±0.01, shown by 
the dashed line. Each point gives the mean amplitude and flux for 
a single observation. The scatter is commensurate with the mean 
uncertainty in fractional amplitude, which is shown by the error 
cross in the upper right. 

the uneven sampling periodicities at 1 d and 96 min due 
to the RXTE observation schedule and orbit. There are 
no peaks at the 2 hr binary orbital period, indicating 
that the pulse profile is independent of orbital phase. 

The resulting noise powers are around 2 decades higher 
at long periods («3 d or longer) than at short periods 
for 1998, and even more for 2005. The 2002 outburst 
spectra exhibit less profile noise at long timescales, but 
still are somewhat red. Poisson statistics produce an 
uncolored lower limit on noise. This white noise domi- 
nates at timescales shorter than the orbital period, ex- 
cept in the case of the particularly noisy fundamental of 
2005. The spectra of the intrinsic profile noise (i.e., the 
spectra after subtracting off the Poisson contribution) 
roughly followed a power law noise spectrum, which we 
parametrized as Pk{f) cx; J-tpln. The best- fit values of 
7PLN, listed in Table [21 varied from roughly 0.4 to 1. 

3.4. Fractional amplitudes of the harmonics 

In our time-domain discussion of the pulse profiles, an 
apparent trend is the tendency of the pulses to become 
narrower, more asymmetric, or doubly peaked — gener- 
ally speaking, to become less sinusoidal — as the out- 
burst's flux decreases. In the frequency domain, the re- 
lation is striking: the fractional amplitude of the 802 Hz 
second harmonic, r2, strongly anticorrelates with the 
background-subtracted 2-25 keV flux, /x, as shown in 
Figure [5] This power-law dependency has a slope of 
—0.50 ± 0.01. The agreement with the data is excellent 
for such a simple model, giving a reduced statistic of 
x1 = 1.15 with 1816 degrees of freedom. It spans two and 
a half decades in luminosity and includes every detected 
harmonic amplitude from all four outbursts. 

In terms of the pulse profile, the second harmonic con- 
tributes in two ways. If its peak is 45° out of phase with 
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TABLE 2 

Noise properties of the outbursts 



Fundamental Second harmonic 

(o-J)^/2 (Tl,int 7PLN {^lY^^ f^Z.int 7PLN 

1998 Apr 0.007 0.014 0.96 ± 0.09 0.021 0.017 0.43 ±0.12 

2000 Feb 0.023 0.052 — 0.022 0.039 — 

2002 Oct 0.012 0.016 0.67 ± 0.09 0.023 0.027 0.51 ± 0.09 

2005 Jun 0.013 0.061 0.85 ± 0.07 0.019 0.024 0.77 ± 0.05 

Note. — All phases are in cycles (i.e., fractions of the 2.5 ms spin period). 
^(T^)^/-^ gives the mean contribution of Poisson noise; <yk;int is the amplitude of 
pulse profile variability in excess of the Poisson noise; and 7pln is the slope of 
the power law best fit to the spectrum of . We did not attempt to estimate 

power law noise slopes for the 2000 outburst because of its low-quality data. 



the peak of the fundamental, it wiU produce an asym- 
metric pulse profile (e.g., profile 6A in Fig. [3]). If it is 
in phase, a narrower primary pulse with a small second 
peak will result (as in profile 6B) . If the components are 
90° out of phase, the profile will be profile 6B flipped, 
but we never observed such a configuration. 

To further understand the influence of flux on the pulse 
profile, we decomposed the second harmonic's fractional 
amplitude into its asymmetric and double-peaked com- 
ponents, 

^2 , asym = r2 | siu 47r'(/' | and ( 1 1 a) 

r2,dp = r2 |cos47rV'| , (lib) 

where tp is the phase offset between the peaks of the two 
harmonics: tp = (02 + ^4'2) — (0i + A(/)i). The resulting 
plots have substantially more scatter than Figure [3] due 
to the uncertainty of ip, which is considerable, particu- 
larly at low fluxes. However, they both roughly conform 

— 1/2 

to the r2 oc /x power law. We conclude that the de- 
crease in flux increases the asymmetry of the pulses and 
the presence of secondary pulses in approximately equal 
measure. 

In contrast, the fractional amplitude of the fundamen- 
tal behaves unpredictably. During the slow-decay stage 
of 1998, it is unvarying and strong, at a constant 5.5% 
rms. During this stage of 2002, it is weaker (4%) and 
somewhat variable; during 2005, it is weaker still and 
erratically changing by up to a full percent between ob- 
servations. Its behavior is more consistent in the tail. 
In all outbursts, the fractional amplitude of the funda- 
mental varies widely, usually (but not always) having its 
maxima around the peaks of the flares and its minima 
during the fading portion of the flares. 

For the most part, a pulse profile model only includ- 
ing the fundamental and second harmonic adequately 
describes the folded profiles. However, folding long 
stretches of data does sometimes result in the detection 
of a third harmonic with fractional amplitudes ranging 
up to «0.25% rms. We do not rehably detect any higher 
harmonics. 

3.5. Upper limits on the subharmonics 

With some assumptions, we can strongly constrain the 
presence of subharmonics and half-integral harmonics. 
The most straightforward approach is to fold all the ob- 
servations using multiples of the best-fit frequency mod- 
els from each outburst. The amplitude of the result- 
ing profile will give an upper limit. The resulting 95% 



TABLE 3 
Upper limits on subharmonics and 
half-integral harmonics 



Harmonic 


Upper limif 


(% rms) 


Factor 


Hz'' 


A 


B 


C 


1/4 


100.2 


0.017 


0.019 


0.52 


1/2 


200.5 


0.022 


0.024 


0.45 


3/2 


601.5 


0.018 


0.021 


0.43 


5/2 


1002.4 


0.026 


0.024 


0.42 



These background-corrected upper lim- 
its are quoted at the 95% confidence level. 
These limits result from combining all the 
observations (column A), combining only 
bright observations (B), and not combin- 
ing any observations (C). See the text for 
more details.'' Frequencies listed here are 
approximate. The upper limits were ob- 
tained using exact multiples of the best-fit 
constant-i' models. 

confidence upper limits are listed in column A of Ta- 
ble [S] However, this approach is only statistically valid if 
the uncorrected fractional amplitude (i.e., the fractional 
amplitude relative to the source counts and the back- 
ground) is constant. Clearly this assumption is false. 
Aside from the varying proportion of source photons, the 
background-subtracted fractional amplitudes of the fun- 
damental and the second harmonic fluctuate throughout 
the outburst, spanning nearly an order of magnitude in 
the tails of the burst. There is no reason to believe that 
a subharmonic would not fluctuate similarly. 

Column B of Table [3] takes the more moderate ap- 
proach of only folding together observations during which 
the 2-25 keV flux exceeds 5x 10^^" erg cm^^ s^^, thereby 
only including the main body of the outbursts. Back- 
ground photons are thus a much smaller contribution, 
and the fractional amplitudes of the observed two har- 
monics were relatively stable during these times. Nev- 
ertheless, we still are folding enough photons to obtain 
very stringent upper limits: in the case of the 200 Hz 
subharmonic, we get a 95% confldence upper limit of 
0.024% rms. We feel that these numbers are our most 
reliable, not making unreasonable assumptions about the 
fractional amplitude fluctuations. 

For completeness, we also include the most conserva- 
tive upper limits, which make no assumptions whatso- 
ever about the fractional amplitudes of the subharmon- 
ics. For instance, it would be possible in principle for 
the subharmonic to be present only during a single ob- 
servation and zero-amplitude everywhere else. To con- 
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TABLE 4 

Best-fit constant frequencies, and their, v upper limits 
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Fig. 6. — Constant-frequency measurements of the 
SAX J1808.4— 3658 outbursts, showing the spin down of the star. 
The frequencies are relative to vq = 400.97521000 Hz. The error 
bars are estimated using Monte Carlo simulations of phase resid- 
uals with the same noise properties as the actual outburst; they 
do not account for the uncertainty in the source position. The 
x's mark what the frequencies would be if the fit source position 
differed from the actual position by 2 tr along the ecliptic plane 
in the direction of increasing RA. The same position error in the 
decreasing RA direction would move the frequency points by an 
equal amount in the opposite sense. 

strain the resulting upper limits at least somewhat, we 
again only used observations during which the source was 
brighter than 5 x lO^^*' erg cm~^ s~^ and that had at 
least 10^ counts. These single-observation limits are tab- 
ulated in column C. 

The stringent upper limits of column B provide the 
best evidence yet that the the spin frequency of the star is 
indeed 401 Hz. If the star was spinning at 200.5 Hz, with 
two antipodal hot spots each emitting pulses to produce 
the observed frequency, a 200.5 Hz subharmonic would 
almost certainly be present. 

3.6. Spin frequency measurements and constraints 

We initially performed the simplest possible fits to 
the phase residuals of each outburst: constant-frequency 
models. We did not include the data at the very begin- 
ning of the 2002 and 2005 outbursts, where pulse profile 
changes during the rise and peak obscure any variations 
in the phase. We also excluded the residuals during 
2002's mid-outburst pulse profile change, but included 
the residuals of the fundamental on both sides of the 
shift by using different profile templates before and after 
it. 

The resulting frequency measurements are shown in 
Figure inland summarized in Tabled These data clearly 
indicate that the source is spinning down. The probabil- 
ity that the actual spin frequency is constant or increas- 
ing is less than 10~^ given the uncertainty estimates. 
These uncertainties do assume that our optical position 
is exact, but the position error is excluded because its 
effects are highly correlated; for instance, the 1998 and 
2002 outbursts are six months apart on the calender, so 
a position offset would produce equal and opposite fre- 





Data included 










(MJD) 


(mHz) 


(10- 


14 Hz s-l) 


1998 Apr 


50914.8 - 50936.9 


0.371 ±0.018 


(- 


-7.5,7.3) 


2000 Feb 


51564.0 - 51601.9 


0.254 ±0.012 


(- 


-1.1,4.2) 


2002 Oct 


52565.0 - 52602. 8=: 


0.221 ±0.006 


(- 


-1.3,2.5) 


2005 Jun 


53529.6 - 53581.5 


0.195 ±0.016 


(- 


-0.5,2.4) 



^ The frequencies are relative to uq = 400.975210 Hz.'' 95% con- 
fidence intervals from the Monte Carlo simulations.'^ Excluding 
MJD 52575.7-52577.7. 



quency displacements for the measurements from these 
outbursts. There is no position that would provide a 
statistically feasible constant or increasing frequency. 

The linear fit through the measured frequencies is not 
particularly good: its statistic is 9.7 with 2 degrees 
of freedom, yielding a probability of about 1% that the 
frequencies are drawn from a linear progression. Once 
again, changing the source position does not significantly 
change the result or improve the fit, and changes in the 
position by more than the 1 a uncertainty along the eclip- 
tic substantially worsen the linear fit. To estimate the 
uncertainties of the linear slope in light of this poor fit, 
we rescaled the measurement errors such that reduced 
statistic would be unity. The resulting first-order spin 
derivative is j> = (— 5.6± 2.0) x 10~^^ Hz s~^. The large 
1 a uncertainty reflects the uncertainty in the slope of the 
frequency change, not in the observation that the source 
is spinning down. The probability that the frequency is 
not decreasing is less than 10~^, as mentioned above, a 
confidence of better than 6 a. 

Fitting second-order frequency models established that 
!> is consistent with zero during all the outbursts. These 
measurements are particularly sensitive to pulse profile 
variations, so care must be taken to not overfit such fea- 
tures. We again exclude the initial observations of the 
1998, 2002, and 2005 outbursts, because the pulse pro- 
file changes would induce large non-zero i/ measurements 
that most likely do not reflect the spin of the under- 
lying neutron star. Using tempo to find the best-fit 
!>'s and applying Monte Carlos to estimate their uncer- 
tainties, we arrived at the 95% confidence intervals of 
Table [H Excluding the 1998 outburst, which had the 
shortest span of timing data and thus the most poorly 
constrained j>, these 95% confidence upper limits were all 
of order \v\ < 2.5 x IQ-^^ Hz s^^. 

The uncertainties in the measurement of the frequency 
preclude phase connection between outbursts. During 
the 920 d gap between the 2002 and 2005 outbursts, the 
6 nHz frequency uncertainty from 2002 would accumulate 
to a phase uncertainty of 0.5 cycles; the 2 x 10^^^ Hz s~^ 
uncertainty in the long-term spin down would contribute 
0.6 cycles. Worse, these estimates are best-case scenar- 
ios, since they assume that the spin down is constant. 

During the 1998 and 2002 outbursts, we observed an 
abrupt change in the slope of the phase residuals at the 
end of the main outburst. We modeled these apparent in- 
stantaneous changes of frequency by including frequency 
ghtchcs in our tempo fits. (While the tempo glitch 
models are useful in describing the data, we do not be- 
lieve that we observed actual sudden changes in the spin 
frequency of the star, a point discussed in detail in ^14.31 ) 
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Fig. 7. — Comparison of the 1998 and 2002 glitch-like events. 
The phase plots show the phase residuals relative to a constant- 
frequency model for the fundamental (black points) and the har- 
monic (grey points), binned such that there is one point per ob- 
servation. The black lines indicate the best timing models fit by 
TEMPO. The middle plot shows the 1998 and 2002 light curves 
for comparison. The 1998 light curve has been vertically offset by 
1 erg cm~^ s~^ for clarity. The data are displayed such that the 
apparent changes in frequency are aligned. Notice that this align- 
ment also has the effect of closely matching up the light curves. 

Figure [7] shows the phase residuals of these outbursts 
and their best-fit ghtch models. These models only em- 
ploy an instantaneous change in frequency; including a 
phase jump or introducing a v after the events did not 
significantly improve the fits."'^^ 

In both outbursts, these apparent frequency changes 
coincide with the sudden drop in flux that marks the 
transition from the slow-decay stage to the flaring tail 
stage. At the same time, the fractional amplitudes of the 
fundamental and harmonic increase, and, in the case of 
2002, the pulse profile change occurs. (This pulse profile 
change, discussed earlier in ii3.21 is apparent in Fig. [7|as 
the rapid advance of the fundamental phase.) If we view 
the phase residuals with respect to the pre-transition fre- 
quencies, as is the case in Figure [71 the residuals following 
the transition skew upward, indicating progressively in- 
creasing lags. This effect is more pronounced in 1998, but 
its coverage is far better in 2002. If we were to interpret 
these changes in slope as abrupt spin frequency changes, 
they would represent drops of 0.21 /iHz and 0.03 /iHz 

During the 2002 outburst, the absence of a phase jump refers 
only to the second harmonic, which we believe is a better tracer 
of the neutron star spi n during this period of time (in agreement 
with the conclusions of lBurderi et al]|2006l) . 



for 1998 and 2002, respectively. (Again, we consider this 
scenario unlikely; see §4.31 ) If we instead interpret them 
as the motion of a radiating spot, the drift rates would be 
6.5° d~^ and 1.0° d~^, retrograde. The total observed 
shifts between the start of the flaring tail and the loss 
of the signal are substantial: 0.15 cycles (54°) in 1998 
and 0.06 cycles (22°) in 2002. The data are not good 
enough to distinguish whether these drifts are continu- 
ous. For instance, it is possible that the hot spot made 
a retrograde jump every time there was a flare. 

We did not observe the main body of the 2000 out- 
burst, so we cannot measure whether the apparent fre- 
quency decreased when it entered the flaring tail stage. 
But if it did, and if the decrease in the apparent frequency 
was of similar magnitude to that observed in 1998 and 
2002, then including the main body of the 2000 outburst 
would raise the overall frequency of the outburst some- 
what. This correction might put it in line with the other 
frequency measurements in Figure [51 reducing the large 
statistic of the constant- fit. Therefore we cannot 
conclude that the change in the observed frequency from 
one outburst to the next is incompatible with a linear 
progression. 

During the 2005 outburst, the substantial pulse pro- 
file noise during the tail prevented us from measuring 
a change in apparent frequency. The uncertainty in 
the measurement of the frequency during the tail was 
0.03 /zHz, as estimated using Monte Carlo simulations of 
the profile noise, and the phase residuals jumped by as 
much as 0.1 cycles from one observation to the next. If 
there was a smaller drift, as seen during the 2002 out- 
burst, we would not necessarily detect it. 

3.7. Evolution of the binary orbit 

We fit the orbital parameters separately for each out- 
burst. Table [5l lists the results. As expected, the val- 
ues of axSini and Porb were consistent among the out- 
bursts. The fit parameters esinw and ecoscj were con- 
sistent with zero. We used them to improve significantly 
on previous upper limits on the eccentricity. 

The measured time of ascending node advanced with 
each outburst, relative to the times expected if the pe- 
riod was constant. Figure [51 shows these Tasc residuals. 
A quadratic provides a good fit (x^ — 1.01 with a sin- 
gle degree of freedom), yielding a constant orbital pe- 
riod derivative of Porb = (3-5 ± 0.2) x 10~^^ s s~^ and 
a significance of 15.6 cr. In a,n inde pendent analysis of 
the same data, iDi Salvo et al.l (|2007f ) report a consistent 
value for Porb- They derive a smaller uncertainty and 
larger x^, most likely reflecting an underestimate of the 
orbital phase measurement errors. 

Table [6l summarizes all the parameters for the pulse 
timing of SAX J1808.4-3658. 

4. DISCUSSION 

Our analysis of multiple outbursts from 
SAX J1808.4— 3658 allows us to greatly improve 
our understanding of the behavior of this low-mass 
X-ray binary. By comparing the observed frequency 
from each outburst, we can see the long-term spin 
down, which is too small to be detectable from a 
single outburst. Comparison of the pulse profiles from 
each outburst lead us to conclude that we are seeing 
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TABLE 5 

Binary parameter measurements from each outburst 



Outburst 


(s) 


ax sin i 
(light-ms) 


Tasc 

(MJD, TDB) 


esinoj 
(10-6) 


e cos LU 
(10-6) 


1998 Apr 
2000 Feb 
2002 Oct 
2005 Jun 


7249.1553(18) 

7249.1565(6) 
7249.1547(24) 


62.8080(46) 

62.8147(31) 
62.8282(109) 


50921.7584194(12) 
51591.8019861(40) 
52570.0186514(9) 
53524.9944192(32) 


-60 ± 64 

8 ±57 
-173 ±83 


-86 ± 64 

41 ± 57 
53 ± 83 



Note. — We excluded the 2000 outburst when calculating everything but Tasc because 
its data were noisy and sparse. 
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Fig. 8. — Measurement of an orbital period derivative. The 
points show the observed times of ascending node, relative to the 
expected times for a constant period. The Tasc of each outburst 
comes progressively later, indicating a period derivative of (3.5 ± 
0.2) X 10-12 s s-l. 



TABLE 6 

Combined timing parameters for SAX ,11808.4-3658 



Orbital period, Porb (s) ^ 

Orbital period derivative, Porb (10-^^ s s^^) 
Projected semimajor axis, ax sin i (light-ms) 
Time of ascending node, Tasc (MJD, TDB) 
Eccentricity, e (95% confidence upper limit) 
Spin frequency, u (Hz) ^ 

Spin frequency derivative, u (lO-^^ Hz s" -'^) 



7249.156961(14) 
3.48(23) 
62.8132(24) 
52499.9602477(10) 
< 1.2 X 10-** 
400.975210240(11) 
-5.6(2.0) 



Porb S'lid !^ a-rG specified for the time Tasc. 

characteristic, repeated profile changes as the outbursts 
progress, rather than a purely random noise process. 
Finally, fitting of the orbital parameters over the seven 
years of observation provides a greatly improved orbital 
ephemeris. 

4.1. Long-term spin down 

By observing the mean spin frequency of each outburst, 
we found that SAX J1808.4— 3658 is spinning down at a 
rate of ;> = (-5.6 ± 2.0) x lO^^^ Hz s-^. This spin 
down results in a loss of rotational energy at a rate of 
E = Att^Ivv — % X 10'^^ erg s~^, assuming a canonical 
value of / = 10"'^ g cm^ for the neutron star (NS) moment 



of inertia. 

Most of this spin down occurs during X-ray quiescence; 
accretion torques during the outbursts play a minimal 
role. Over the seven years of our observations, the mean 
outburst frequency decreases by 1^2005 — '^iggs = — 0.18± 
0.02 /iHz. Let us suppose that this frequency change 
happens only during the X-ray outbursts (which have 
a duty cycle of <5%). Since the outburst light curves 
are quite similar, it is reasonable to presume that each 
would contribute roughly the same frequency shift, thus 
splitting this frequency change into three equal steps. If 
the spin down is due to a constant z>outburst that acts 
during the «20 d of each outburst, -'^'^ then 



-0.18 AiHz _,4 

I'outburst ~ „n J ^ 5 X 10 Hz S 

3 X 20 d 



-1 



(12) 



By contrast, we were able to set stringent (95% confi- 
dence) upper limits of |z>| < 2.5 x 10^^* Hz s-^ during 
the outbursts (Table H]). We conclude that the spin down 
is dominated by torques exerted during X-ray quiescence. 

We will thus consider three possible sources of torque 
during quiescence: magnetic dipole radiation, the expul- 
sion of matter by the magnetic field (i.e., the propeller 
effect), and gravitational radiation. In general, we as- 
sume that all three mechanisms contribute addititively 
to the observed spin down of SAX J1808.4-3658, 

A^obs = Afdipolc + Afprop + iVgr. (13) 

We discuss each below. 

4.1.1. Magnetic dipole torque 

A spinning dipolar magnetic field will produce a sig- 
nificant spin down during quiescence for the 10* G 
field strengths expected for a millisecond pulsar. Rel- 
ativistic f orce-free MHD rn odels of pulsar magneto- 
spheres by ISpitkovsicvl (|2006f) give a torque of A^dipoio = 
~^'^{2ttv/cY{1 + sin^ a), where /i is the magnetic dipole 
moment and a is the angle between the magnetic and 
rotation al poles. Pulse profile rnodelin g of the 1998 out- 
burst bv lPoutanen fc Gierhnskil (|2003f ) suggests that the 
magnetic hot spot is not far from the rotational pole, sep- 
arated by an angle of 5-20°. While other effects might 
also contribute to the spin down, the rotating magnetic 
field will always be present and provides an upper limit 

-'^ In reality, i> would almost certainly not be constant during as 
the accretion rate varies, but for argument's sake we make the most 
conservative assumptions possible. A varying i> would req uire that 
it be sometimes greater than the value from equation (|12^ , making 
it even less plausible that it would escape detection. 
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on the dipole moment: 



H<0.77 X 10^6 {l + sin^a) ^''^ 

1/2 



10*^ g cm^ 

— V 



401 Hz 

1/2 



-3/2 



5.6 X 10^16 Hz s- 



G cm-" 



(14) 



For CL — 15°, this upper limit on the dipole is 0.75 x 
10^^ G cm'^, yielding a field strength of roughly B — 
1.5 X 10^ G at the magnetic poles. We emphasize that 
this upper limit on the magnetic field is for a purely dipo- 
lar field. The presence of higher-order multipoles would 
require a stronger field at the NS surface to produce the 
observed v. This field estimate is consistent with the 
limits implied by accretion physics (see ii4.5p . 

If magnetic dipole torque is a significant contribu- 
tor to the spin down of SAX J1808.4-3658, then the 
source may behave like a rotation-powered pulsar dur- 
ing quiescence, producing radio pulsations and a par- 
ticle wind. The heating of the companion by a par- 
ticle wind has been invoked as an explanation of why 
the compan ion is significantly br ighter than expected in 
the optical. iBurderi et al.l (|2003f) predicted a dipole mo- 
ment of /i = 5 X 10^^ G cm'^ based on the optical ob- 
servations, somewhat higher than our approximate up- 
per limit on /x, but most likely within the uncertain- 
ties of the model. A similar analysis by iCampana et al.l 
(|2004f ) found the needed irradiation luminosity to be 
-^x = (4!ti) X 10^'^ erg s~^, compatible with the observed 
E — % 10'^^ erg s~^ loss of rotational energy. No radio 
emission has been de tected during quiescence. The up- 
per limits of 0.5 mJy fGaen sler et al.lll999l :[B 



urgav 



et al.l 



12003:!) are not particularly constraining. 

The X-ray luminosities of isolated millisecond pulsars, 
for which magnetic dipole radiation is the primary spin- 
down mechanism, shows a strong correlation with their 
rates of ro tatio nal energy loss . From the t ables com- 
piled in iZavlini (2006 ) and Cameron et all (|2007f l. the 
5-10 keV X-ray luminosity goes as oc E}-^^ with less 
than a quarter decade of scatter. Based on this empiri- 
cal relation, we would expect a quiescent luminosity for 
SAX J1808.4-3658 of 5 x 10^° erg 8"^ However, this 
prediction is a factor of ten lower than the observed qui- 
escent fluxes of 8 x 10^^ erg s~^ and 5 x 10'^^ erg s~^ 
(jCampana et al.l 120021 : iHeinke et all I2007D . suggesting 
other mechanisms for quiescent emission are at work. 

4.1.2. Magnetic 'propeller torque 

The propeller effect offers another possible explanation 
for the observed spin down during quiescence. If the Kep- 
lerian corotation radius [tco = \GM j ^■k^v^^I'^ « 31 km) 
is less than the magnetospheric radius tq, at which point 
the infalling matter couples to the magnetic field, then 
the magnetic field will acc elerate the matter, possibl y 
ejecting it from the system (jlUarionov fc Sunvaevlfl97l . 



The lSpitkovakyj I I2006I ') formula for Af^ipoie differs substantially 
from the classically derived torque due to a rotating dipole in a 
vacuum, A^vac = ^/t'^(27ri//c)^ sin-^ a, especially for small a: for 
a = 15° , the derived limit is approximately one fifth of the vacuum 
value. 



The torque exerted on the neutron star by propeller ejec- 
tion of matter at a rate Moj depends on the details of the 
interaction between the pulsar magnetosphere and the 
accretion disk. However, we can parametrize this torque 
as 



prop ' 



-nAfej(GMro)i/2 
= -n(ro/rco)^/V/ej(GMrco)^/' , (15) 

where the detailed physics determines the dimensionless 
torque n, which is zero for rp = r^a and of order unity 
for To > 1.1 (|Eksi et al.l l2005V 

We can then roughly estimate the rate at which matter 
would need to be ejected from the system during quies- 
cence to account for the observed spin down: 



Afcj < -2.3 X 10 



-12 



^(ro/rc 



-1/2 



lO'^s g cmV \\A M( 



M 



-2/3 



(401 Hz) 



5.6 X 10 



-16 



Hz s- 



Mq yr- 



1/3 



(16) 



As a consistency check, we note that this upper limit 
does not exceed the predicted long-term mass transfer 
rate for the binary, 1 x 10~^^ Mq yr~^, which is driven 
by gr avitational radiation emission due to the binary 
orbit (|Bildsten fc Chakrabartvl I2001D . Indeed, not all 
the mass lost by the donor star will necessarily reach 
the pulsar magnetosphere during quiescence and be pro- 
pelled outward; most of it would queue up in the ac- 
cretion disk and later r each the NS during an outburst. 
iGallowav &: Gumming] (|2006D found that the mass trans- 
fer is roughly conservative, albeit with enough uncer- 
tainty that propeller mass loss as large as the above Moj 
limit is not ruled out. 

Even if propeller spin down provides the dominant qui- 
escent torque, the resulting ejection of matter from the 
system would not greatly affect the binary orbit. The 
timescale for propeller spin down is proportional to the 
timescale for the ejection of mass: Porb/-Porb (x Mej/Mc, 
where « 0.05 is the mass of the companion. 
Applying the above Mej gives Mc/Me] = 20 Gyr, far 
longer than the observed orbital evolution timescale of 
^orb/-forb = 66 Myr. More refined calculations using 
the arguments of iTauris fc van den Heuvel ()200d) yield 
a propeller timescale of 6 Gyr, still far too large. Glearly 
there are other contributions to the orbital evolution; we 
discuss some in fJH 



4.1.3. Gravitational radiation torque 

A variety of mechanisms have been proposed in 
which rapidly rotating neutron stars can develop mass 
quadrupoles that give rise to gravitational radiation 
from the neutron star itself. These mechanis ms include 
r-mo de instabilities (Wagoner 1984; Andersson et al.l 
|1999( ). acc retion-induced variations in the density of the 
NS crust (|Bildstenl 119981 : lUshomirskv et al.l |200i QD, dis- 
tortio n of the NS due to toroidal magnetic fields (jGutled 
120021 ). and magnetically confined m ountains at the mag- 
netic poles ()Melatos fc Pavnj|2005D . The loss of angu- 
lar momentum due to gravitational radiation has been 
suggested as a mechanism to explain the absence of ob- 
served pulsars with spin frequencies faster than ^730 Hz 
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(|Chakrabartv et al.|[200l IChakrabartvll2005f) and makes 
millisecond pulsars a target for interferometric gravita- 
tional wave detectors. 

The mass quadrupole moment of the star, Q, deter- 
mines the torque produced by gravitational radiation: 
Ng,- = —^GQ'^{2-Ki>/c)^. For our measured this sets 
an upper limit of 




oi' Q ^ 10^* /. The strain amplitude of the resulting 
gravitational waves, averaged over all NS orientations, is 
he = llbGv^Q / dc^ (Brady et al. 1998), giving a charac- 
teristic strain at Earth of /ic = 6 x 10~^^. This strain 
is undetectable by current or planned gravitational wave 
experiments. For Advanced LIGO, with a strain sensi- 
tivity of ^3 x 10"^^ Hz-^/^ in the 100-400 Hz range 
(iFritschell 120031 ). even a search using an accurate phase 
model would require years of integration time. Note that 
the dependence of N^-^ on the v is very strong, so it is 
quite possible that gravitational wave emission produces 
larger spin downs in faster («700 Hz) rotators. 

4.2. Pulse profile variability 

The evolution of the pulse profile is clearly not purely 
stochastic. With multiple outbursts, we are able to note 
for the first time that the pulse profile seems to take 
on similar shapes at similar times in the outbursts, as 
illustrated in Figure [3l These characteristic changes in 
the pulse profiles suggest that the emitting regions of 
the NS are changing shape and position as the outbursts 
progress. The consistency of these changes, along with 
the consistency of the outburst light curves, suggests that 
as the accretion disk empties onto the star, the geometry 
of the disk, the accretion funnels, and the resulting hot 
spots evolve in a similar manner for each outburst. 

The most striking example of ordered pulse-profile evo- 
lution is the strong relationship between the harmonic 
content and luminosity: r2 oc L"^!"^ . Given the complex- 
ity of the system, its abidance by such a simple model 
is quite surprising. SAX J1808.4— 3658 is not alone in 
this behavior. At least two other millisecond pulsars, 
IGR J00291-H5934 and XTE J1807-294, exhibit similar 
inverse correlations between the amplitude of their sec- 
ond harmonics and luminosity (Hartman et al. 2007, in 
prep.). 

One possible explanation is recession of the accretion 
disk as the accretion rate drops, revealing the star's pre- 
viously occulted second hot spot. For rapidly rotating 
pulsars, partial occultation of the star by the accretion 
disk will be common. Assuming a mass of 1.4 Mq, the co- 
rotation radius of SAX J1808.4-3658 is Tco = 31 km « 
3i?, where R is the NS radi us. Following the sta ndard 
pulsar accretion model fe.g.. lGhosh &: Lambl ll979). the 
inner edge of the accretion disk will be at roughly the 
Alfven radius: ro « = (2G'A/)-i/'^Af-2/7^4/7 rpj^jg 
truncation radius must be at rg < r^a for infalling matter 
to reach the NS surface. There are clear problems with 
the application of this model, which was developed for 
higher-field pulsars with ro ^ R: the width of the transi- 
tion region in which the magnetic field becomes dominant 



is on the same order as its distance to the star, muddling 
the definition of a truncation radius. Nevertheless, this 
simple model is still qualitatively instructive. 

Since r^o ~ 3i?, neutron stars in systems with incli- 
nations i > 70° will always be partially occulted during 
outburst. During the outbursts of SAX J1808.4-3658, 
the peak fluxes at which the pulses are most sinusoidal 
are roughly a factor of 10 greater than the low fluxes 
at which the harmonics are more prevalent (cf. Fig. [5]). 
As a result, the Alfven radius will increase by a factor 
of ?'A,taii/''A,poak ~ 10^^^ « 2 as the source dims. Be- 
cause the maximum Alfven radius is «3i? during accre- 
tion, the radius during the peak of the outbursts must 
be rA,pcak ^ %R- At this separation, the star will be 
partially occulted if i > 45°. Thus the degree of occulta- 
tion will depend on M for a wide range of inclinations. 
For 45° < i < 70°, the disk wiU partially occult the 
NS above some critical M. For i > 70°, the NS will 
always be partially occulted, with the degree of occulta- 
tion increasing as M increase s. Pulse proflle modeling by 
IPoutanen &: Gierhnskil (|2003f ) suggests that the system is 
at an inclination of i > 65°. 

The observations that show clearly double-peaked 
pulse proflles happen exclusively in the final, flaring tail 
stage of the outbursts, typically during the fading por- 
tion of a flare. In view of this model, one could imagine 
that the accretion disk is most recessed as the flares fade. 
One difficulty with this model is that the increased r2 
observed at low luminosities is not solely due to the ap- 
pearance of doubly peaked pulse profiles; many profiles 
in this regime show single pulses, but with substantially 
greater asymmetry than typically seen at higher lumi- 
nosities. 

Another possible cause is the expansion of the hot 
spots during high accretion due to diffusive effect s . Sim - 
ulations of accretion flows by iRomanova et al.l ()2004l ) 
demonstrate that as the fluence increases, the cross- 
sections of t he ac cretion funnels grow. Modeling by 
iMuno et al.l (|2002f l establishes that the harmonic con- 
tent of the pulsations decreases as the size of the hot 
spot increases. 

4.3. Motion of the hot spot 

During the 1998, 2002, and 2005 outbursts, we ob- 
served clear trends in the phase residuals that sug- 
gest that the emitting regions do not remain at a flxed 
longitude. During 2002 and 2005, an abrupt phase 
change in the fundamental at the end of the main 
body of the outburst produces an advance of the pulse 
peak that corresponds to a shift of the hot spot by 
«50° eastward. These shifts are simultaneous with 
and occur on the same 3-4 d timescale as the sud- 
den drops in luminosity at the end of the main out- 
bursts. During 1998 and 2002, the phase residuals of 
both harmonics begin gradually increasing during the 
flaring tails of the outburst, corresponding to a west- 
ward drift of the hot spots. Motion of the hot spot 
has also been suggested to explain ph ase residuals in 
GX 1+4 and RX J08 12.4-3114 (iGallowav et al.l [200l 
and XTE J1814-338 (jPapitto et al.ll2007f) . 

For a more natural description, we adopt the Earth-based 
convention of longitude: earlier pulse arrivals = prograde hot spot 
motion = eastward shift, and vice versa. 
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These trends in the phase residuals almost certainly 
represent motion of the observed hot spot rather than 
frequency glitches. Glitches are rapid changes in the spin 
frequency of the NS due to imperfect coupling between 
the crust a nd more rapidly rotati ng, supcrfluidic lower 
layers (e.g.. lAnderson fc Itohlll975h . This interaction oc- 
curs well below the accretion layer, and it would not be 
expected to coincide with or have the same timescale as 
rapid changes in the accretion rate. 

When discussing the motion of the hot spots, the lon- 
gitudes of the magnetic poles provide natural meridians 
from which to measure phase. Since their movement 
would require the realignment of currents in the core and 
crust, the magnetic poles remain at fixed positions for 
timescales far longer than the outbursts. The suppres- 
sion of regions of t he field due to accret ion also occurs 
on long timescales (jCumming et al.|[2001[ ). 

For high-field pulsars, the magnetospheric radius is far 
from the star, and the accretion column follows field lines 
that reach the NS surface near the magnetic pole. This 
is not necessarily the case for low-field pulsars. A closer 
accretion disk will intersect more curved field lines, which 
terminate farther from the poles. In the previous section, 
we described how the Alfven radius can move outward 
from roughly 1.5R to 3R as the accretion rate drops. 
As the disk recesses, it will intersect decreasingly curved 
field lines that are rooted closer to the poles, causing the 
hot spots at the bases of the accretion columns to also 
approach the poles. 

This simple picture can explain the observed phase 
shift as the luminosity rapidly drops during the end of 
the 2002 and 2005 outbursts. In both cases, the lumi- 
nosity decreases by about a factor of 4. A change in 
M by this magnitude would cause the Alfven radius to 
move outward by a factor of 1.5 and the inner edge of 
the accretion disk to move outward by a similar amount. 
This change will almost certainly cause material removed 
from the inner edge of the disk to attach to a different 
set of field lines, with the larger radius favoring lines that 
attach closer to the pole. If the hot spot tends t o be to 
the west of the pole, as seen in iRomanova et al.| (|2004[ ) 
for a magnetic pole an angle of a = 30° from the rota- 
tional pole, then the attachment to different field lines 
would produce an eastward drift as observed. That said, 
these MHD simulations appear to have strong, chaotic 
dependencies on their parameters. (For a = 15°, the hot 
spot is south of the magnetic pole; for 30°, west; and for 
45°, north!) More work is needed to better model these 
observations. 

This scenario does not explain why the shift of the 
pulse peak would solely be expressed by a change in the 
fundamental; during these episodes in 2002 and 2005 the 
phase of the harmonic remains relatively constant. How- 
ever, a movement of the hotspot toward the magnetic 
pole would most likely change the shape of the hotspot, 
possibly in a way that would preserve the phase of the 
harmonic. 

The slow drifts seen during the tails of the 1998 and 
2002 outbursts are also difficult to explain. The flares 
during the tail cause the luminosity to change in a peri- 
odic manner, so we cannot expect a monotonic motion 
of the accretion disk's inner edge. The net drift during 
the tail of the 2002 outburst is of the same magnitude as 



the rapid phase shift that happens right before the tail 
begins, suggesting the drift may be a relaxation of the 
accretion column back to its original location. 

4.4. Comparison with previous spin 
frequency measurements 

There have been a number of previous re- 
ports of short-term i) measurements made during 
outbursts of several accreting mi llisecond pul- 
sars including XTE J0929 -314 (iGallowav et all 



200l. SAX J180 8.4-3658 (iMorean et all 120031: 
Burderi et all l200l . XTE J1 751-305 (iMarkwardtl 
2003D. ICR J0 0291-H5934 (iFalanga et al.l 120051: 



Burd eri et aLll2007t ). and XTE J1814-334 (jPapitto et all 
2007). Some of the reported i> values have been surpris- 
ingly large given the estimates of M during the outbursts, 
possibly violating a basic prediction of magnetic disk 
accretion theory: that accretion torques cannot exceed 

the characteristic torque iVchar — M{GMrcoY^'^ exerted 
by ac creting Keplerian mat erial at the corotation radius 
(e.g.. lGhosh fc Lambl[l979l) . 

In the particular case of SAX J1808.4— 3658, spin 
derivatives as large as a few times 10~^^ Hz s~ ^ near 
the outburst peak were reported (|Morgan et al.l l2003l : 
iBurderi et al.l 12006} ) . corresponding to accretion torques 
exceeding A^char for this source. However, these studies 
calcula ted pulse phase resid uals using only a single har- 
monic; [Morgan et al.l (|2003D reported v det e ctions using 
only the fundamental, while IBurderi et al.l ([2006') mea- 
suring the phase from the second-harmonic alone after 
noting the sudden phase shift of the fundamental in the 
middle of the 2002 outburst. Our results in §3.6 indicate 
that both of these approaches are likely to be contam- 
inated by pulse shape changes, at least in the case of 
SAX J1808.4-3658. 

Figure [9] illustrates this point. Taking the phases of 
the harmonic components as direct spin measurements 
can produce large values of u during the peak of the 
2002 outburst. Fitting only the fundamental's phase 
residuals during the first 10 d of the outburst, we find 
V = (-1.2 ± 0.4) X 10-13 Hz s-^ On the other hand, 
using only the second harmonic for the same interval, we 
find v = (5.3 lb 0.1) X 10-13 Hz s-^, in good agreement 
with the lBurderi et al.] ()2006f ) measurement. Because the 
pulse shape is changing rapidly during this part of the 
outburst, the pulse arrival times cannot be accurately de- 
termined. We therefore cannot reliably use this part of 
the outburst to measure the spin of the NS. Note that if 
we exclude this region of large pulse shape variability, the 
remaining phase residuals are consistent with a constant 
spin frequency over the outburst interval (§3.6). From an 
examination of all the outbursts of SAX J1808.4— 3658 
(excluding regions of large pulse shape variability), our 
work sets an upper limit of \v\ < 2.5 x lO-^"' Hz s-^. 

We thus conclude that the past measurements of short- 
term v in SAX J1808.4— 3658 are unreliable. The analy- 
sis technique we described in §2.3 can mitigate the effects 
of pulse shape variability to some extent, but attempts 
to measure v in accreting pulsars must properly account 
for these variability effects, and in some instances these 
effects may prevent such measurements. The v mea- 
surements reported in other accreting millisecond pulsars 
must all be reevaluated in this light; all the apparent 
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Fig. 9. — Fitting a frequency model using only the fundamental 
(black) or the harmonic (grey) produces non-zero v measurements 
during the peak of the 2002 outburst. The data points are the 512 s 
phase residuals relative to the best constant-frequency model for 
the 2002 outburst. The solid lines give the best constant-;> models, 
fit solely to the fundamental or the second harmonic. The dashed 
line shows the constant-frequency model derived using both, com- 
bined via equation JOj; this fit did not use the points prior to MJD 
52565. 



violations of the iV < A'char limit predicted by theory 
may be owing to spurious measurements caused by pulse 
shape variability. However, at least some accreting mil- 
lisecond pulsars are observed to have relatively stable 
pulse shapes, indicating that accurate short-term v mea- 
surements are possible and that previous measurement 
of these sources should be reliable. 

4.5. Constraints on the magnetic field 

We showed in §4.1.1 that the condition iVdipoie < ^obs 
implies that the magnetic dipole moment fi < 0.8 x 
10^^ G cm'^. This limit is consistent with the range for 
H implied by the observation of a ccretion-powered pulsa- 
tions throughout the outbursts (jPsaltis fc Chakrabartvl 
|1999( ). At low accretion rates, the field cannot be so 
strong that it centrifugally inhibits matter from reach- 
ing the NS; during times of high accretion, it must 
be strong enough to truncate the disk above the stel- 
lar surface in order for there to be pulsations. The 
dimmest observation in which we observed pulsations 
was in 1998, with a flux in the 2-25 keV band of 
1.5 X 10~^^ erg cm~^ s~^; the brightest was at the peak 
of the 2002 outburst, 2.62 x lO'^ erg cm'^ s^^. These 
fluxes, along with an improved estimate of the Edding- 
ton luminosity fro m observations of photosphe ric radius 
expansion bursts (jGallowav fc Cummingj [2006[ ) . give us 
new limits on the range of accretion rates at which pulsa- 
tions have been detected, relative to the Eddington rate 



Me: M„ 



1.8 X 10-4 Me and M„ 



0.03 Me. 



(We have made the usual assumption that L (x M.) 
These limits allow us to update the range for fi derived in 



iPsaltis fc Chakrabartvl (|1999D . equations (11) and {12):'^° 

0.2 X 10^** G cm^ < ^ < 6 X 10^^ G cm^ . (18) 

Taken together with the A'dipoic limit, wc obtain a fairly 
narrow allowed range for the magnetic dipole moment. 



0.2 X 10^^ Gem-' <fi< 0.8 x 10^^ G cm' 



(19) 



which corresponds to a surface dipole magnetic fleld 
strength of (0.4-1.5) x 10® G. This field is relatively weak: 
the magnetic fields implied by the A ustrilia Telescope 
Natio nal Facility Pulsar Catalog^^ (jManchester et aH 
|2005( ) for millisecond pulsars range from 1.1 x 10® G to 
14 X 10® G. 

4.6. Constraints on accretion torques 

Even though we did not detect an accretion-induced ;> 
during the outbursts, our new upper limits on |;>| provide 
far stronger constraints on the accretion physics of low-i? 
systems such as SAX J1808.4— 3658 than previous mea- 
surement s. Following the earlier a nalysi s of the 1998 out- 
burst bv IPsaltis fc Chakrabartvl (|1999( l. the lower limit 
on the spin frequency derivative predicted by accretion 
torque theory during an outburst with an average accre- 



tion rate of Af, 



-M 



0.01 Me is 



!> > 2 X 10" 



M 



1.4 M 







1045 g 
1/2 



R 



10 km 



3/2 



M 



avg 



0.01 Me 



Hz s" 



(20) 



where 77 is a dimensionless parameter encapsulat- 
ing the disk-magne tosphere interaction. (Refer to 
iGhosh fc Lambl[l979l for a discussion of the physics that 
goes into this parameter.) 77 is strongly dependent on 
the nicigneto spheric radius. For tq ~ t^co^ 

the NS will be 

in spin equilibrium with the accreted matter and rj will 
be small. From the confidence intervals in Table SI 
the probability that we would have missed detecting the 
resulting 2 x 10^^"* jjg s^^ spin up is 0.15%, suggesting 
that rj < I and the source is near spin equilibrium during 
the outbursts. 

4.7. Discussion of the increasing Porb 

Our seven year baseline for timing analysis provides 
the most precise measurements yet of the orbital period 
of SAX J1808.4-3658. We find that the orbital period is 
increasing at a rate Porb = 3.5(2) x 10"^^ s s"^. This Porb 
lie s somewhat ou t side th e 90% confidence upper limit set 
bv iPapitto et al.l (|2005f) using the 1998—2002 outbursts, 
most likely owing to the more limited baseline available 
in that analysis. 

It is interesting to compare our measurement with 
theoretical expectations. For orbital periods <3 hr. 

In derivin g this range for fi, we make t he sa me conservative 
assumptions as lPsaltis fc Chakrabartvl 1)19991 ): the lGhosh fc LambI 
Ill991h boundary layer parameter ranges on 0.1 < 7s(M) < 1; 
the NS mass is 1.4 Mq < M < 2.3 Mq; and the NS radius is 
10 km < R < 15 km. 

http: //www. atnf . csiro . au/research/pulsar/psrcat/ 
Pulsars associated with clusters were excluded to minimize the 
impact of line-of-site accelerations. Field strengths were approxi- 
mated using equation II14I I 
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mass transfer is LMXBs is driven by angular momen- 
tum losses due to g ravitational radiation from the binary 
(jKraft et al.lll962( ). since magnetic br aking torques are 
thought to be ineffective in t his regime (jRappaport et al.1 
[lOM iSpruit fc Ritteii flOS^l ). For SAX J1808.4-3658, 
the M predicted by this mechanism is consistent 
with observationally inferred long-term average value 
of M = 1 X 10"" Mq yr-i fBild sten fc Chakrabartvl 
1200 It) . For conservative mass transfer from a degener- 
ate (brown dwarf) donor, this predicts orbital expan- 
sion on a time scale PgrbZ-Pp rb = 3 Gyr (see, e.g., 
iTauris fc van den Heuvell 120061 ). By contrast, our mea- 
sured value of Poib/-Poib — 66 Myr is an order of magni- 
tude more rapid. 

The origin of the anomalously large Porb in 
SAX J1808.4— 3658 is unclear, although we note that 
unexpectedly large Porb values have also been ob- 
served in several other L MXBs including 4 U 1820—30 
(Ivan der Klis et all [1993 ^. EXO 0748-676 (|Wolff et al.l 
|200^, and 4U 1822-371 (HcUicr ct al. 1991). As pointed 
out by iChakrabartv fc Morgan! (|19980 . the binary pa- 
rameters of SAX J1808.4— 3658 are very similar to those 
of the so-called "black widow" millisecond radio pulsars, 
all of which arc ablating their low-mass companions (see, 
e.g., Fruchter ct al. 1990). If SAX J1808.4-3658 does 
indeed turn on as a radio pulsar during X-ra y quies- 
cence (jBurderi et al.l[2003t ICampana et al.l[2004t see also 
§4.1.1), it may be a black widow system as well, consis- 
tent with its very low donor mass. As such, it is interest- 
ing to note that a large and variable Porbi both positive 



and negative, has been measured in two black widow pu l- 
sars (jArzoumanian et al.lll994l : iDoroshenko et al]l2001h . 

Although mass loss from the companion through 
an ablated wind would tend to increase Porb, 
the mass loss rate required to explain the ob- 
served Arb in SAX J1808.4-3 658 is -lO-^M© yr"! 
(jTauris fc van den Heuvell l2006( l: this is unphysically 
large given our measured pulsar spindown rate (§4.1), 
which sets the pulsar luminosity available for irradiat- 
ing the companion. This explanation for Poib is also 
inadequate in the black widow pulsars, where the or- 
bital period variability is quasi-cyclic on a 10 yr tim e 
scale (jArzou mania n et al.lll994l : IDoroshenko et al]l2001h . 
In those systems, it has been suggested that tidal dis- 
sipation and magnetic activity in the companion is re- 
sponsible for the orbital variability, requiring that the 
companion is at least partiall y non-degenerate, convec- 
tive, and magnetically active (lAr zoumanian ct al.* |1994l : 
lApplegate fc Shahaml[l99 4l: IDoroshen ko et al. 200i[ ). If 
this mechanism is active in SAX J1808.4— 3658, we would 
expect quasi-cyclic variability of Porb to to reveal itself 
over the next few years. 
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APPENDIX 

A. IMPROVED OPTICAL POSITION FOR SAX J1808. 4-3658 

An accurate source position is essential for high-precision pulsar timing. An incorrect position results in errors during 
the barycentcring of X-ray arrival times, producing frequency offsets due to improperly corrected Doppler shifts (see, 
e.g., Manchester fc Peters 1972). SAX J1808.4-3658 hes only /3 = -13.6° below the ecliptic plane, so any errors 
during barycentcring will be particularly pronounced. For example, a position error of e = 0'.'2 parallel to the plane of 
the ecliptic produces frequency and frequency derivative offsets relative to « 401 Hz of 

Ai/ = i/Qe (a© cos/3/c) {2n/P^) cost = 40 cost nHz (Al) 

Ai> = -iyoe{a:^cosP/c) {2Tr/ P^^^f sinr = -8 x 10"^^ sin r Hz s"^ . (A2) 

Here r — 2Tit/PQ parametrizes the Earth's orbit, with time t equal to zero when the Earth is closest to the source. 
These offsets are comparable with the expected timing uncertainties. Each outburst gives a baseline of about 2 x 10^ s 
over which we can typically measure pulse arrival times with an accuracy of better than 25 /^s, or 1 x 10~^ cycles, 
producing ^5 nHz frequency uncertainties. By similar logic, we should be sensitive to v^s as small as ~3 x 10^^^ Hz s~^. 
In practice, the pulse shape noise observed in SAX J1808.4— 3658 makes the actual uncertainties somewhat greater 
than these back-of-the-envelope values, increasing the v uncertainty by a factor of ^2 and the v uncertainty by a factor 
of ~10, but the frequency uncertainty is still substantially less than the offsets due to a 0'.'2 position error. 

We observed the field of SAX J1808.4— 3658 with the Raymond and Beverly Sackler Magellan Instant Gamera 
(MagIG) on the 6.5-m Baade (Magellan I) telescope on the night of 2001 June 13, using the r' filter. The seeing was 
0'.'5. Figure [To] shows the results. After standard reduction, involving bias-subtraction and flatfielding, we attempted to 
register the field to the International Goordinate Reference System (IGRS). We examine d three astrometr i c catalogs for 
this purpose: the Hubble Space Telescope Guide Sta r Catalog (GSC, which was used bv lGiles et al.lll999 !: 'Laske r et al. 
[1990), the USNO-Bl.O survey (|Monet et al.ll2003l) . and the Two-Micron All-Sky Survey (2MASSnSkrutskic eLaL 
I2OO6'). We selected stars from all three catalogs that were not saturated or blended on our image, and fit using the 
IRAF task ccmap for the position offset, rotation, and plate-scale. We found that we could obtain the best astrometry 
with 2MASS: with USNO and GSG, many stars that had consistent positions between 2MASS and our image had 
deviations of more than 0'.'2, the overall scatter was larger, and there were fewer stars. With 2MASS we fit using 70 stars 
across the 2' MagIC frame. With position residuals of 0'.'08 in each coordinate, we obtained a combined uncertainty of 
0'.'08/\/70 = O'.'Ol. Therefore, our astrometric uncertainty is dominated by the R:i0'.'15 position uncertainty of 2MASS. 

To verify our position, we checked for stars on the MagIC image from the Second US Naval Observatory CCD 
Astrograph Catalog (^UCAC2: IZacharias et al.ll2004f ). These are highly accurate positions (individual uncertainties of 
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Fig. 10. — A 30" portion of our r'-band Magellan image. The counterpart of SAX J1808.4— 3658 is indicated by the tick marks: it is 
the north-west object of the close pair at the center. We also indicate three 2MASS stars that we used for astrometry with the circles. 
The changing grayscale levels across the image reflects poor correction for the four- amplifier readout of MagIC but does not affect our 
astrometry. 



20-40 mas) for relatively bright («15 mag) stars taken with a CCD at a current epoch (1996-1998) and with proper 
motions. We found three unsaturated UCAC2 stars on our image: 16259696, 16259777, and 16259680. We measured 
their positions on our image and compared the positions derived from the 2MASS solution to those from UCAC2, 
updated to epoch 2001.45. We found no net shift, and the offsets are less than 0'.'16 in all cases. (We note that the stars 
are toward the edge of the image, where residual image distortions may be present, in contrast to SAX J1808.4— 3658 
which is at the center of the image). Therefore we believe that our solution using 2MASS is indeed accurate to our 
stated uncertainty of 0'.'15. 

We then measured the position of SAX J1808.4— 3658 on the image and transformed the position to the ICRS. The 
position that w e find is: R . A. = 18'^08"^27!62, Decl. = -36°58'43'.'3, equinox J2000.0, with uncertainty 0'.'15. This 
is 1'.'5 from the iGiles et al.l ()199S) pos ition, twice their quoted 0'.'8 uncertainty. But with many more reference stars 
of higher quality over a smaller field ([Giles et al.lll999l used 5 GSC stars over a 4' field), and CCD data taken at a 
more recent epoch (1998 for 2MASS, vs. 1987-1996 for GSC^^ and 1981 for USNO), this new position should be more 
accurate. 

B. DERIVATION OF PHASE UNCERTAINTIES 

Derivation of the uncertainties of the phase residuals, as given in equation (|4]), follows from our definition of the 
phases, 



Ak exp {2mk(j)k) = 2 Xj exp {2mjk/n) , 
where we have divided our phases into n bins, each containing Xj photons. Inverting to solve for 0^, 



(Bl) 



2nik 



jk - In 4^ 



(B2) 



where we define the constants Ejk = exp (2TTijk/n) for the sake of brevity. 

For relatively low fractional amplitudes (certainly the case throughout this paper), each phase bin will contain 
approximately the same number of photons: Xj ~ Np^/n, with variances (crxj)^ « N-ph/n due to Poisson counting 
statistics. These add in quadrature to give the variance in 0^: 



o(Pk \ , .1 



1 



E. 



A^ph 
n 



27rfc X)"'=i Xj'Ej>k J 

The USNO does not recommend GSC for current use: see |http : //ad.usno .navy ■mll/star/star_cats_rec . shtml#gsc2 .2\ 



(B3) 
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Summing the exponentials, we have 



'jk 



From the definition of Aj. in equation (|Bip . ^j-^jk 



^Ak- Substituting these in, we reach our estimate of the phase uncertainty: o-fc = ^^2Nph/2^TkAk 
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